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1 Introduction 

Let A be the generator of a one-parameter semigroup T t acting in a Hilbert space 
TC. We discuss the numerical computation of T t f, or equivalently the solution of 
the initial value problem 

fit) = Af(t) (1) 

given /(0) = /. This involves several problems. The first is that the spectral 
mapping property may fail; that is one may have 

Spec(Ti)\{0} {e xt : A G Spec(A)}. 

In particular ||T t || may grow faster than e s * as t increases, where 

s = sup{Re (A) : A G Spec(^4)}. 

This problem is well known and has been studied from many points of view, but 
it remains a difficulty, even if A has discrete spectrum, E3 EH ECU IJ3 HZ] • 

The second problem arises for differential operators, particularly in several space 
dimensions, when the matrix approximations have very high dimensions. Even 
if A has a sparse matrix, T t generally has a full matrix, so storing the matrix 
entries is not feasible. The obvious solution is to find a subspace of relatively small 
dimension which contains most information of interest. One might try to do this 
by taking the linear span of a finite number of eigenvectors, those for which the 
real parts of the eigenvalues are largest. Unfortunately experience shows that for 
many non-self-adjoint operators A, the eigenvectors do not form a basis; indeed 
the norms of the spectral projections often increase exponentially fast according 
to natural orderings of the eigenvalues. This forces one to be very cautious about 
assuming that a spectral expansion of some given / will yield useful results. 

This phenomenon is linked to the appearance of non-trivial pseudospectra. When 
this happens the determination of more than a small number of eigenvalues may 
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become numerically impossible. Even if theorems about the convergence of the 
eigenfunction expansion of a general / G B subject to a resummation method can 
be proved, they have limited use if most of the eigenvalues and eigenvectors cannot 
be determined. 

Many recent papers about pseudospectra have drawn attention to possible instabil- 
ity problems which are not revealed by looking at the spectrum alone, HI EH 
HHII2DI- Our goal in this paper is more positive: we use pseudospectral methods to 
solve the evolution equation above for highly non-self-adjoint operators. The exis- 
tence of a large number of approximate eigenvalues is regarded as a resource rather 
than an embarrassment. We develop an 'approximate spectral expansion' which 
may have little to do with the true eigenvalues and eigenvectors of the operator. 
In spite of this our main result, Theorem may be used to solve the evolution 
equation to a high degree of accuracy. In the examples studied numerically we 
demonstrate that it is far more accurate than the normal spectral expansion. 

Our method is particularly useful if one wishes to solve the initial value problem 
(fT|) for a large number of different choices of the initial data. The approximate 
eigenvalues and eigenvectors only need to be produced once, and the computations 
needed for each choice of the initial data are fairly easy. 

The examples which we consider in this paper are convection-diffusion operators. 
There are arguments in favour of studying the associated semigroups T t = e At in 
L 1 rather than L 2 . Diffusion is a probabilistic phenomenon, and the conservation 
of probability is not easy to study in an L 2 context. It is shown in [21 El M EE] 
that the 'same' semigroup may have different growth properties when studied in 
L 1 or in L 2 . Nevertheless we will focus on the L 2 theory, for the same reason as in 
classical Fourier theory: the theorems are much simpler to state and apply. 



2 The Abstract Setting 

We start with several assumptions. The first is the choice of numbers M, 7 such 
that 

\\T t \\ < Me 7i (2) 

for all t > 0. The second is the existence of a set S equipped with a a-field of 
subsets and a finite measure ds. We assume that we are given measurable families 
of unit vectors u s G Ti and of complex numbers X s parametrized by s G S and 
satisfying 

inf{||-u s — to|| + \\Aw — X s w\\ : w G Dom(v4)} < e. (3) 

Throughout this paper e is a given 'acceptable' error satisfying < s < 1/2. From 
a purely theoretical point of view the assumption 

\\Au s - \ s u s \\ < e (4) 
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for all s G S would be simpler. We prefer (J3J) because it permits simpler expressions 
for the vectors u s in applications. Clearly \ s are approximate eigenvalues of A, up 
to the error e > 0. The assumption @ implies that 

A s G Spec 2e (A) :={z:\\ {zl - A)' 1 1| > (2s)- 1 } 

in the language of pseudospectral theory. 

If A is highly non-self-adjoint, the fact that X s are approximate eigenvalues of A 
does not imply that they are close to the spectrum of A. This allows us to go far 
beyond what is possible by means of conventional spectral methods. In numerical 
applications we will take S to be finite, but the above setting allows a better 
understanding of the general theory. 

We define a bounded, linear 'pseudospectral' transform Q from L l (S) to Ti by 

G<fi= 4>(s)u s ds. 
Js 

We restrict Q to L 2 (S) and note that it is then bounded with ||£?|| < l^ 1 / 2 , where 
\S\ is the measure of S. The adjoint operator Q* : Ti — > L 2 (S) is given by 

(07)00 = </,«.> 

and B = G*g : L 2 (S) -> L 2 (S) is given by 



W)(s)= / b(s,t)<f>(t) dt 



s 



where 



b(s,t) = (u t ,u s ). 



Since b is a bounded measurable function, B is a Hilbert-Schmidt operator on 
L 2 (S). 

It is immediate from the definitions that B(p = if and only if Qcf) = 0. We assume 
throughout the paper that B is invertible, a matter which needs to be confirmed 
in any application. 

The following theorem describes how best to approximate / G Ti. by expressions 
of the form Qcj) where <fi G L 2 (S). We will frequently refer to the algebraic sum 
M. = £ + where C is the range of Q in TC. This is a dense linear subspace of 
TC. If S is finite, as in all numerical applications, then £ is closed and M. — Ti. 

Theorem 1 If P is the orthogonal projection on Ti with range L, then 

Pf = QB~ x Q*f 

for all f G M.. For such f the quantity \\Q<f> — f\\, where <p £ L 2 {S), is minimized 
by (f) = B v Q*f. We also write 4> = Q\f, as in Matlab. 
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Proof If / G £ x then Q* f = 0, so QB~ l Q*f = 0. If / = Q<f> then 

QB~ x G*f = Q{B- l Q*Q)$ = G<P = f. 
This proves the first statement. If / = Qifj + g where if) G L 2 (S) and g G C then 

ll/-M a = IIW-^)ir + IMI 2 . 

This is clearly minimized for <f> = if) and, under our standing hypothesis that Q is 
one-one, this is the unique minimum. We also have 

B~ l g*f = B- l g*{Qif) + g)=ifj. 

The above method of approximation should be contrasted with the following alter- 
native. Suppose that A has a complete set of eigenvectors u n , n = 1, 2, ... and that 
m* are corresponding eigenvectors of A*, so that the two sets form a biorthogonal 
system in the sense that {u ni u* m ) = S m>n . The standard spectral expansion with 
respect to this system is 

/ = lim Q N f (5) 

where 

N 

QnJ = 51 (/,<)«„. (6) 

n=l 

If the identity (JHJ) holds for all / G 7i one says that u n form a basis in TC. Un- 
fortunately this is rarely true for highly non-self-adjoint operators. Indeed \\Qn\\ 
frequently diverges at an exponential rate in applications. One might modify the 
above formula by assuming Cesaro or Abel summability, but convergence would 
still have to be verified and is not always true. 

One the other hand if the set {m„}^ ! L 1 is complete we always have 

/ = hm P N f 

where is the orthogonal projection of 7i onto lin{w n : 1 < n < N}, and 
this is indeed the optimal approximation sequence to /. Putting S = {1,...,N} 
the theorem above enables one to compute Pn- The main disadvantage of the 
projections P/v is that they do not commute with A. 

Returning to the general context at the start of this section, we use the operators 
defined above to solve the evolution equation approximately. We start by obtaining 
a bound on the real parts of approximate eigenvalues. 

Lemma 2 If ||w|| = 1 and 

\\u - w\\ + \\Aw - Xw\\ < e 

then Re (A) < 7 + 2Mb. 
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Proof We first observe that < e < 1/2 implies 1/2 < ||w|| < 3/2. Putting 
\x = Re (A), the identity 



implies that 



{T t _ s e Xs w} = T t _ s e Xs (\w - Aw) 



\\T t w-e xt w\\ = || fr t _ s e As (Aw- Aw) ds\ 
Jo 

< C Me 7(t - s)+fiS e ds 
Jo 



= Me- 



/i-7 



It follows that 



e M 72 < 2Me 7 * + Me 



e^ - e 7 * 



/i-7 

If /x > 7 then letting t — > +oo, we deduce that 

2Me 
1 < 



/i-7 

which is equivalent to the statement of the lemma. 

It follows immediately from the lemma that if Re (A) > 7 and we put A = 7 + 
ilm(A) then |A - A| < 2Me. Therefore 

\\u - w\\ + \\Aw -Xw\\ < e(3M + 1). 

In the rest of the paper we assume that these changes in the approximate eigen- 
values have been made, so that Re (A s ) < 7 for all s e S, and that e has been 
increased correspondingly. 

The main theorem of this paper is best formulated in terms of certain approximat- 
ing semigroups R t . 

Theorem 3 Let e xt be the multiplication operator on L 2 (S) defined by 

(e xt <P)(s) = e x ^(s). 
and define R t on M. for t > by 

R t = g e xt B- 1 g* (7) 

Then R = P, R t {C L ) = and R t {£) C C for all t>0. We also have 

RtRuf — Rt+uf 

for all t,u>0 and f e M. 
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Proof If / e C L then Q*f = so R t f = 0. If / = Q§ where G L 2 (S) then 

RJ = gePB-^g^ = g(e xt <f)) e C. 

Finally, if / = g<f> then 

R t R u f = i^(e A » = g(e xt e Xu <P) = R t+U f. 

Theorem 4 Suppose that S is finite and 7i = L 2 (X, dx), and rewrite u s (x) = 
u(x,s). Then 

(RJ)(x)= [ K t (x,y)f{y) dy 
Jx 

for all f EH, where 

K t {x,y) = $>(x, s)e Xst (B- l ) s , r u(y,r) 

r,s 

Proof Since S is finite, £ is a finite-dimensional subspace of H, and C + C 1 = H. 
We deduce that Rt has domain 7i. The formulae of the theorem are the result of 
rewriting (J7J) in integral operator form. 

One might conjecture that the integral kernel of Rt is uniformly close to that of T t 
under suitable conditions, but we do not have any such results. 

The following is our main theorem. It is only numerically efficient if e > and 
5 = ||/ — Pf\\ are both small. We discuss this further in the next section. 

Theorem 5 If f <E M then 

\\T t f-Rtf\\ < \\f-Pf\\Me^ + e{l + M + Mt)\\g\f\\ x ^ (8) 

for allt > 0. 

Proof If we put <fi = g\f = B~ l g* f then the estimate can be rewritten in the 
form 

\\T t f-g<p t \\ < ||/-e0||Me^+ £ (l + M + Mt)||0|| ie ^ (9) 
where <p t = e xt (p. We follow the argument of Lemma El up to 

\\T t w - e xt w\\ < f Me 7(i - s)+ " s £ ds 
Jo 

< [ Me 7t e ds 
Jo 

= eMte 11 . 

Hence 

||T t M-e A 'M|| < \\e xt {w-u)\\ + \\T t {u - w)\\ + \\T t w - e xt w\\ 
= e(l + M + Mt)e^. 
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Applying this to each u s in the expansion 

G<j>= <j>(s)u s ds 
Js 

yields 

\\T t G<t> ~ Gfa\\ <s(l + M + Mt)||0|| ie ^. 
The theorem follows by combining this with the bound 

\\T t f-T t G4>\\<Me«\\f-G4>\\. 

The above theorem is only useful as long as the right hand side of (JHJ) is much 
smaller than ||-Rt/||. Since 

WW = \\GU < UWi^ 

where 

/i = sup{Re (A s ) : s G S}, 

the estimates are only useful for a short time if fj, <C 7. The point here is that 
\i may be substantially larger than sup{Re (z) : z G Spec(A)}, so pseudospectral 
methods may be correspondingly more accurate than spectral methods. 



3 Numerical implementation 



In numerical applications we take S to be a finite set, possibly containing fewer 
than a hundred points. This implies that M. = TC. The main task is the choice 
of the vectors u s G TC. Once this has been done, there are three possible methods 
of computing = B~ 1 G*f given / G TC. The vectors u s determine G*, and also 
the operator B via the kernel b(s,t). One might compute B~ l and then apply the 
above formula to obtain 0. Since the operator B" 1 is highly singular it is better 
to evaluate B^ip for ip = G*f without computing B~ l \ Matlab uses the command 
B\ip for this purpose. One may finally avoid any reference to G* or B, by using 
Matlab to compute <fi = G\f directly. Since G is a rectangular matrix, Matlab 
actually finds the 'solution' with least squares error. We tried all three methods, 
and found, as expected, that the third is by far the most accurate. Once <fi has 
been determined we do not use Theorem El as stated, but the reformulation in ©. 

The choice of a suitably small e > is made before starting the computation. On 
the other hand the verification that ||/ — G$\\ is small is done on a posterior basis. 
Since and G have to be computed in any case, this poses no problems. 

There are two obvious ways of choosing the unit vectors u s . One may use one's 
physical intuition, as in the examples of this paper, to select certain vectors, and 
then show that they satisfy the fundamental inequality (jH)) for a suitably small e > 
0. This method has been used successfully in the semiclassical context, |H1 IH El E2j ■ 
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The second method depends upon numerical, pseudospectral calculations, and will 
be described in more detail in a later publication. The first stage is the replacement 
of the differential operator A by a sparse matrix approximation, possibly in a space 
of very high dimension. This may involve finite element methods or wavelets, and 
is not the focus of this article. There is now a well-developed technology for 
calculating pseudospectra, and it may be applied to very large sparse matrices. 
Given e > 0, we next have to choose a finite set of numbers from the set Spec £ (v4). 
If t > is known, there is no need to consider points in A G Spec e (yl) for which 
e A * is extremely small, because the contributions of the corresponding terms of Q<pt 
will be negligible. This applies in particular to any eigenvalues of A whose real 
parts are much less than 7. For each \ s we finally choose a unit vector u s for which 
\\Au s — \u s \\ < e. 

In some cases it might be advisable to chose several vectors u s corresponding to 
each A s , providing each vector with a different label s. The choice depends upon 
how many eigenvalues of order e 2 the operator 

D s = (XJ-Ay(X s I-A) 

possesses. For rotationally invariant problems in dimension two, for example, one 
would treat each angular momentum sector independently, and include that pa- 
rameter in the labelling of S. The e-pseudospectra for different sectors may well 
overlap. 



4 A Pure Convection Operator 

The theory above has applications to convection-diffusion operators, but the sim- 
plest example is given by the pure convection operator 

(Af)(x) = f(x) 

acting in L 2 (0, a) subject to the boundary condition f(a) = 0. This is the generator 
of the one-parameter semigroup T t given by 




f(x + t) if x + t < a 
otherwise. 



Since T t = for all t > a, 7 can take any value in the estimate (j2J- Nevertheless, 
since we are interested primarily in the case of large a, we take 7 = and M — 1. 
The fact that Spec (A) = implies that one cannot hope to use spectral expansions 
to evaluate T t , but pseudospectral expansions are still possible. Since this example 
is exactly soluble, we only analyze it by our method in order to understand how 
well the method works. We will see in the next section that the pseudospectral 
expansion of this operator is an asymptotic form of the corresponding expansion 
for a simple convection-diffusion operator. 
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The following constructions depend upon the choice of positive constants c and a. 
In many cases an appropriate value of c may be found in the range 5 < c < 10. 
The value c = leads to a Fourier series expansion, which is not appropriate for 
this problem. One could put a — 1, but for asymptotic theorems it might be more 
appropriate to make it proportional to a and/or inversely proportional to c. Let 
v : [0, 1] — > [0, 1] be the function 



v\x) 



1 if0<x<a — a 

(a — x)/a if a — a < x < a. 



(Many other choices would be equally suitable, for example v(x) = 1 — e^ x a ^ a .) 
Given s G Z we define u s G L 2 (0, a) by 

uJx) = ke~ cx/a+27risx/a 



where 

a 



2c 

This choice implies the identity = 1. We also see that ah 2 /2c — > 1 at an 
exponential rate as c increases. If we define w s G Dom(y4) by 

w s (x) = u s (x)v(x) 

then 

\\u s ~w s \\ 2 = k 2 [ a e~ 2cx/a (l - v(x)) 2 dx 
Jo 

nk 2 

~ 2c 

If we put A s = —c/a + 2nis/a then 

\\Aws- \ s w s \\ 2 = k 2 [ a e~ 2cx/a v'(x) 2 dx 

Jo 

ak 2 



< 



-2c(l-a/a) 



2ca 2 

This indicates that the bound (jSJ) holds with e = 0(e _c(1_Qv/a )) as c — > oo. 
Having chosen a sufficiently large N > 0, we then put 

S = { s G Z : -N < s < N}. (10) 

The integral kernel of B = Q*Q is 

b(s,t) = (u t ,u s ) 

= k 2 / e _2cx / a + 27ri (* _ ' s ) :r / a qI-j; 







l-e" 2c 



2c/a — 2 r Ki{t — s)/a 
1 



1 — izi{t — s)/c 
9 



if c is sufficiently large. 

If / e L 2 (0, a) and g = Q* f then 

g (s) = k f a f(x)e- cx/a - 2msx/a dx 



and 

N 

(Pf)(x) = ke~ cx/a (B' 1 g)(s)e 2msx/a . (11) 

s=-N 

If Pf is approximately equal to / then we have shown that 

N 

(T t f)(x) ~ ke~ cx/a Y, (B- 1 g)(s)e 27Tisx/a+i - c/a+2nisx/a)t (12) 

s=-N 

for t > 0. In numerical implementations one actually uses the equivalent formula 

N 

(T t f)(x) ~ fee - "/ £ 0( s ) e 2*W«H-(-c/cH-2*wz/a)t ^3) 
s=-iV 

where = (?\/, in the notation of Mat lab. 

Let us compare this with what one gets by using ordinary Fourier series, by making 
the choices c = and N = oo in the above formulae. We then have k = a -1 / 2 and 

u B (x) = a - 1/2 e 2nisx/a 

for all s G Z. We also have 

g( s ) = a - 1/2 [ a f(x)e- 2msx/a dx 
Jo 

so g is the sequence of Fourier coefficients of /, assuming periodic boundary con- 
ditions. Since 

^ | otherwise, 
B is the identity operator on / 2 (Z), and (fTTjl is replaced by 

oo 

f{x) = a' 1 ' 2 g(s)e 2nisx/a , 

s=—oo 

while (JT2|l is replaced by 

oo 

(f t f)(x) = a~ 1/2 Y g(s)e 2nisx/a+27Tist/a 

s=— oo 

= f{x + t), 

subject to periodic boundary conditions on [0, a]. The use of Fourier series therefore 
solves a different problem from that in which we are interested. 
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We implemented the above ideas numerically for two choices of the initial function. 
We put a = 20 and divided each unit interval into 50 equally spaced points, so 
that functions on [0, a] are approximated by sequences with 1000 terms. We chose 
the initial function to be 

/(x)=2e- 10 ^- 5 ) 2 -e^- 5 ) 2 / 10 . 

We defined f t to be the right-hand side of (fT2"j) and computed 

P=\\f-fo\\oo q= \\T t f -A||oo 

for various values of c, N, putting t = 5. (Similar results are obtained using the 
L 2 norm.) The results are presented in Table 1. Our conclusion from the data is 



c 


N 


P 


1 


5 


30 


0.056 


0.056 


10 


30 


0.049 


0.049 


5 


40 


0.0075 


0.0075 


10 


40 


0.0063 


0.0063 


3 


50 


0.0040 


0.050 


5 


50 


0.00063 


0.0067 


10 


50 


0.00051 


0.00051 






Table 1 





that the errors depend more upon the number of terms 2N + 1 in the expansion 
than upon the value of c. However, in the best case, N = 50, we see that c needs 
to be substantially bigger than 5 for accurate results. 

We also considered the initial function / = 1, for which T t f is the characteristic 
function of (0, a — t). We made the same choices a = 20, N = 50, c = 10 and t — 5 
as above. This case is highly singular, since neither / nor T t f are close to being 
in the domain of A. Although the computed values of f t are close to 1 for x < 15 
and close to for x > 15, there is a Gibbs-type phenomenon near x = 15, the 
maximum value of f t being about 1.21. As expected, the maximum is unchanged 
for N = 100. 

The example of this section may be described in terms of a global approximating 
semigroup, to be contrasted with the local approximating semigroups of TheoremEl 
We introduce the operator 

(A c f)(x) = f(x) 

acting in L 2 (0,a) subject to the boundary conditions f(a) = e~ c /(0). If c > 
this is the generator of a one-parameter semigroup T Cjt acting on L 2 (0,a). One 
sees immediately that u s , X s are the eigenvectors and eigenvalues respectively of 
A c . Section 2 provides estimates of how closely solutions of /'(£) = Af(t) are 
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approximated by solutions of f'(t) = A c f(t) which involve only a finite number of 
eigenvectors of A c . However, the right-hand side of (|12|) is not simply a spectral 
expansion of T Cjt . The similarity between T t and T Cjt explains why we should expect 
steadily better approximations as c increases, provided the computations remain 
feasible. 



5 A Convection-Diffusion Operator 

The difference between the L 1 and L 2 behaviour of semigroups is well illustrated 
by the pure convection operator 

(Af)(x) = -2xf'(x) 

which generates the semigroup 

(TJ)(x) = f(e- 2t x). 

This is a positivity preserving contraction semigroup on Cq(R), but on L 2 (TVj we 
have ||Ti|| = e* for all t > 0. The semigroup has the same behaviour when acting 
on L 2 (—a, a), and if b > is large enough one would expect similar behaviour for 
the convection-diffusion operator 

(Af)(x) = b- 1 f"(x)-2xf'(x) 

acting in L 2 (—a,a) subject to Dirichlet boundary conditions at ±a. 

We consider the somewhat simpler operator 

(Af)(x) = b- 1 f"(x) + f(x) 

acting in L 2 (0, a) in more detail. The first term produces a diffusion effect while 
the second cause a drift to the left at speed 1. If we impose Dirichlet boundary 
conditions then T t = e At is a positivity preserving contraction semigroup on L p (0, a) 
for all 1 < p < oo. If b is large then the norm of T t remains close to 1 for t up to 
about a and then decreases rapidly towards 0. We put M = 1 and 7 = in our 
theorems. 

The following results are well-known, |TT^ fTj . The eigenvectors and eigenvalues of 
A are given by e n (x) = k n e~ bx ^ 2 sm(irnx / a) and A n = —b/4 — ir 2 n 2 /ba 2 respectively 
for n = 1,2 We have ||e„|| = 1 for all n if 



e -bx 



1 

4 Jo 

27r 2 n 2 (l - e- 
b(b 2 a 2 + Air 2 n 2 ) 



^ninx/a ^—irinx/a 



ba \ 



2 

dx 



We see that 



2 



,2 b 3 a 
27r z n z 
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as 6 — > oo for fixed n, a. The spectrum of A is asymptotically empty as b — > oo for 
fixed a. It converges to (— oo, —6/4] as a — > oo, but this is not the spectrum of A 
considered either in P 2 (0, oo) or in L 2 (R). The normalized eigenvectors of A* are 
e* = k n e b ( x ~ a ^ 2 sm(irnx/a). 

Lemma 6 The two sets of eigenvectors {e n } and {e^} satisfy (e n , ej^) = if 
m n. The corresponding spectral projections P n of A satisfy 

\\P e / 

b a a a 

as 6 — ► oo /or eac/i n, a. 

Proof The first statement can be verified directly, but it is a consequence of the 
fact that the two sets are eigenvectors of A and A* respectively. A direct calculation 
shows that 

k 2 a 

(e n ,e* n ) = 2 J a/2 - (15) 
The second statement now follows by substituting (jT3| and (|T3jl into 



LP, 



it 



All of the above facts suggest that one should not use spectral expansions for large 
6. 

In order to test this we computed Pjv/ as defined by Theorem [T] with C = 
lin{ei, eTv} and Qjv/ as defined by ©. We chose a = 20 and discretized using 10 
points per unit interval, so that [0, a] was replaced by a set of 201 points, including 
the endpoints. We took the function / to be 



/(*) 



-{x-a/2) 2 



Table 2 shows the sizes of p — \\f — P/v/|| and q — \\f — QnIW for a range of 
choices of N when 6 = 2.5 and when 6 = 5.0. We see that both methods have 
comparable accuracy for iV = 100. However, the method using P/v attains this 
accuracy far more rapidly as iV increases than the pure spectral method using 
Qn. As 6 increases the convergence of both methods deteriorates, and for 6 = 7.5 
neither method gives useful results for any value of N up to 100. 

Our goal in the remainder of this section is to demonstrate that pseudospectral 
expansions are useful for much larger values of 6. For any choice of 6 the pseu- 
dospectra behave in an interesting way as a increases. For every z inside the 
parabola a G R — > —b~ l a 2 + io~ one has 

lim \\(zl - v4) _1 || = +oo 

a— >oo 

and one can construct approximate eigenfunctions for all such z by the following 
method. Given 5 satisfying < 5 < 1/2 and a G R, we put 

Ua (x) = k ^- b / 2+bS+l ^ X - e (-*/2-M-*)*) (16) 
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h — 


2.5 




h — 


5.0 


N 

i V 


V 




n 
H 


N 

i V 


V 




rt 
H 


i n 


QO v 1(1" 


-l 


1 s y in 3 

l.O A 1U 


i n 

1U 


7 ^ v i n - 

( .O A 1U 


-i 


q q v in 9 

O.O A ±U 


on 
zO 


4.0 X 1U 


-2 


1.3 x 10 3 


on 
zO 


1.1 x 1U 


-i 


1.4 x 10 7 


30 


1.7 x 10" 


-3 


5.7 x 10 1 


30 


4.8 x 10" 


-3 


3.3 x 10 7 


40 


1.8 x 10" 


-5 


1.4 x lO" 1 


40 


1.4 x 10- 


-4 


2.8 x 10 5 


50 


5.3 x 10- 


-8 


2.0 x lO" 3 


50 


1.9 x 10' 


-4 


8.3 x 10 2 


60 


4.1 x 10" 


11 


2.9 x 10~ 6 


60 


3.4 x 10- 


-5 


1.7 x 10° 


70 


1.7 x 10" 


11 


4.4 x 10~ 10 


70 


4.8 x 10- 


-4 


6.0 x 10~ 5 


80 


9.0 x 10- 


12 


1.4 x 10~ 10 


80 


2.4 x 10- 


-5 


5.9 x 10~ 5 


90 


1.2 x 10" 


11 


2.0 x 10~ 10 


90 


1.4 x 10- 


-5 


9.5 x 10~ 5 


100 


9.3 x 10- 


12 


3.4 x 10~ 10 


100 


1.0 x 10- 


-1 


2.1 x 10- 4 



Table 2 



where k = k(b, 5, a, a) is given by 



k 



-2 



a e (-b/2+b5+ia)x _ Q {-b/2-b5-ia)x 2 







dx. 



Clearly \\u a \\ = 1. We make 5 depend upon a according to the formula 

5 = 1/2- c/(ab) 

where < c < ab/2. (As before one might choose c in the range 5 < c < 10.) 
This choice of 5 ensures that k~ 2 ~ a(l — e~ 2c )/2c, |w CT (a)| ~ ke~ c , u a (0) = and 
u' a (0) ~ k(b + 2ia) as a — > oo. We also put 

w a (x) = u a {x)v(x) 

where v(x) — 1 — e( x ~ a )/ a ^ and a is a constant such as a = 1. Finally we put 

\i a = r 1 (-6/2 + 65 + iaf + (-6/2 + 65 + i<r) 
= (65 2 - 6/4) - 6-V 2 + 2i5a 
= -b~ l a 2 +ia- c/a + c 2 /(a 2 b) - 2iac/(ab) 

-> -6 _1 CT 2 + Z(T 

as ci — > oo. 

There is no reason to expect that taking a large value of 6 should cause problems. 
Indeed, as 6 — > oo, the functions u a defined by ()lb|) converge to the corresponding 
functions u s defined for the pure convection operator of Section 4. In both cases the 
size of the constant c controls the degree of accuracy of the fundamental estimate 
©• As we have seen before, this has to be weighed against the increased difficulty 
of performing the computations for large c. 
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Theorem 7 Under the above conditions there exists a constant K a ^ a such that 

IK - w a \\ + \\Aw a - fi a w a \\ < K aM a 1/2 | e2c _ 1 1 
for large enough a > 0. 
Proof We have 

\\u a -W a \\ 2 = k 2 /"° e (-6/2+W+*)* _ e (-b/2-bS-ia)x 2 e 2(x-a)/a ^ 

Jo 

< 4k 2 r e (-b+2b8)x+2( X -a)/a ^ 

JO 

- ^! ^-2c _ e ~2a/a\ 

1/a - c/a " ) 

< 3k 2 ae~ 2c 

for large enough a > 0. 

Since b~ x u" a + u' a — fi a u a and w a G Dom(A), we have 

Aw a — fi a w a = 2b~ 1 u' a v' + b~ x u a v" + u a v'. 

Therefore 

\\Aw a — /io-Wo-H < 2b' 1 \\u' a v'\\ + 6" 1 \\u a v "|| + ll-Uo-f'H. 

Each of the terms on the right-hand side is estimated in the same way as above. 
For example 



I 'Il2 ; 2 -2 

\U a V || — k Qt 



f 

Jo 



c (-b/2+bS+ia)x _ e (-b/2-bS-ia)x 



2 



C 



2(x-a)/a dx 



< 2fc2 fo-2c _ p -2a/a> 

" a 2 (l/a-c/a)V e 



2„ -l„-2c 



for large enough a > 0. Combining all these estimates yields the statement of the 
theorem. 

For general values of a G R the functions u a do not satisfy any set of linear 
boundary conditions. However, if we put a = 2ns /a where s G Z then there 
exist non-zero constants q such that u a (0) = 0, u' a (0) — C\ + c 2 <7, Wo-(a) = c 3 and 
= C4 + c 5 o". Therefore the functions u a all satisfy boundary conditions of the 
form u(0) = and 

C5-u'(0) — C2u'(a) = cqu(o). 

We tested the above ideas numerically. We re-parametrized by means of the sub- 
stitution a = 2ixs/a where s G Z and —N < s < N. We put a = 20, each unit 
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interval in [0, a] being represented by 10 equally spaced points. We put a = 1, 
b = 20 and c = 5. We took the same function / as before, that is 



f( X ) =e -i-a/2)\ 

Table 3 shows the values of p = \\f — -P/v/H f° r various values of N. The dimension 
of the subspace £ is 2N + 1. 



2N + 1 


P 




11 


3.3 x 10" 


-i 


21 


3.4 x 10' 


-2 


31 


1.1 x 10- 


-3 


41 


1.1 x 10' 


-5 


51 


3.4 x 10" 


-8 


61 


2.9 x lO - 


-11 


71 


1.3 x 10' 


-14 



Table 3 



The superiority of this method of expansion over both of the previous ones is 
immediately clear. Further computations show that the pseudo-spectral method 
works just as well for all values of b from 5 to 100 (and probably beyond that). 

We finally computed the approximation f t = Q<j> t to T t f given by the formula (JBJ) 
of Theorem |3J We chose the parameters and initial value of / as above but put 
N = 15; the choice iV = 30 gave the same results up to the accuracy displayed. 
We discovered, as expected, that f t is approximately non-negative; in fact 

-3 x 10" 4 < min{/ t (x) : < x < a} < 

for all < t < 16, at which point we stopped the computation. The shape of 
ft remained approximately gaussian as t increased, with the centre moving to the 
left and the width slowly increasing. The maximum of ft decreases slowly up to 
t ~ 10, when the centre of the peak approaches the origin, after which it decreases 
rapidly. The graphs of /, / 4 , / 8 , /i 2 are plotted in Figure 1. 
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Figure 1. Graphs of /, f$ (dotted) and f^, fyi (solid) 
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The detailed behaviour of the maximum m is presented in Table 4 for c = 10. The 
values are the same for c = 5 up to t = 14 after which they decrease more slowly. 
We compare m with = (1 + 4:t/b)- 1/2 . This is the ' same' constant calculated 
using Fourier transforms when a = oo, i.e. for the semigroup on L 2 (R) when the 
initial function is f(x) = e~ x . The two agree up to t = 10, which is all that one 
could expect. All of the results confirm that the pseudospectral approximation to 
the semigroup is highly reliable for the stated values of a and b, at least for this 
choice of the initial function /. 



t 


m 


moo 





1.0000 


1.0000 


2 


0.8451 


0.8452 


4 


0.7454 


0.7454 


6 


0.6742 


0.6742 


8 


0.6202 


0.6202 


10 


0.5593 


0.5774 


12 


0.1268 


0.5423 


14 


0.0049 


0.5130 


16 


0.0000 


0.4880 



Table 4 



We repeated the calculations leading to Figure 1, but with the initial function 
g(x) = 1 for all x G [0, a]. This is a much more serious test of the method since 
g does not satisfy the boundary conditions even approximately. With N = 15 
and c = 5 we obtained the results shown in Figure 2. One sees that g t is close to 
the characteristic function of [0, a — t], but smoothed out because of the diffusion 
term in A. By contrast with the similar calculation in Section 4, there is no Gibbs 
phenomenon, presumably again because of the diffusion term. For smaller values 
of 6, such as b = 5, the fact that g t (0) = for allt > is much more obvious. 
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1.2 




_0 2 1 1 1 1 1 1 1 1 1 1 > 

' 2 4 6 8 10 12 14 16 18 20 



Figure 2. Graphs of g, g$ (dotted) and g 4 , gi 2 (solid) 
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Table 5 lists the first few eigenvalues A n and approximate eigenvalues fi s of A in 
decreasing order of their real parts, where a = 20, b = 20 and c = 5. The largest 
eigenvalue —5.001 controls the asymptotic decay of the semigroup as t — > oo, but 
it has little influence on the size of ||T t /|| for t = 10. One of the main reasons for 
the accuracy of the pseudospectral expansion is the fact that there are so many 
approximate eigenvalues whose real parts are close to zero. For c = 10 the real 
parts of these /x s decrease from —0.488 to —0.729. 
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